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1. Introduction. Nonlinear dispersive partial differential equations (PDEs) play 
an important role in applications since they appear in many approximations to systems 
in hydrodynamics, nonlinear optics, acoustics, plasma physics, Bose-Einstein conden- 
sates among others. The most prominent members of the class are the celebrated 
Korteweg-de Vries (KdV) equation and the nonlinear Schrodinger (NLS) equation 
and higher dimensional generalizations of these. In addition to the importance of 
these equations in applications, there is also a considerable interest in the mathemat- 
ical properties of their solutions. It is known that nonlinear dispersive PDEs without 
dissipation can have dispersive shock waves [?], i.e., regions of rapid modulated oscil- 
lations in the vicinity of shocks in the solutions to the corresponding dispersionless 
equations for the same initial data. Thus solutions to dispersive PDEs in general will 
not have a strong dispersionless limit as known from solutions to dissipative PDEs as 
the Burgers' equation in the limit of vanishing dissipation. An asymptotic description 
of these dispersive shocks is known for certain integrable PDEs as KdV [?, ?, ?] and 
the NLS equation for certain classes of initial data [?, ?, ?]. For KdV an example is 
shown in Fig. for details see [?] . 

No such description is known for 2 + 1-dimensional PDEs. In addition solutions 
to nonlinear dispersive PDEs can have blowup, i.e., in finite time a loss of regularity of 
the solution with respect to the initial data. It is known for many of the PDEs under 
consideration when blowup can occur, but for the precise mechanism of the blowup 
often not even conjectures exist. 



*We thank P. Matthews, B. Muite, A. Ostermann, T. Schmelzer, who provided example codes, 
and L.N. Trefethen, who interested us in the subject, for helpful discussion and hints. This work has 
been supported by the project FroM-PDE funded by the European Research Council through the 
Advanced Investigator Grant Scheme, the Conseil Regional de Bourgogne via a FABER grant and 
the ANR via the program ANR-09-BLAN-0117-01. 

tinstitut de Mathematiques de Bourgogne, Univcrsitc dc Bourgogne, 9 avenue Alain Savary, 21078 
Dijon Cedex, France (christian.kleinOu-bourgogne.fr) 

flnstitut de Mathematiques de Bourgogne, Univcrsitc de Bourgogne, 9 avenue Alain Savary, 21078 
Dijon Cedex, France (kristelle.roidot9u-bourgogne.fr) 



1 




Fig. 1.1. Numerical solution of KdV for the initial data uo{x) = 1/cosh'^a: (left) and the 
corresponding asymptotic solution (right) for t = 0.35 and e = 10~^ (see [7]). 

In view of the importance of the equations and the open mathematical questions, 
efficient numerical algorithms are needed to enable extensive numerical studies of the 
PDEs. The focus of the present work is to study a 2 + 1-dimensional generalization of 
the KdV equation, the Kadomtsev-Petviashvili (KP) equation, and a 2+ 1-dimensional 
generalization of the NLS equation, the Davey-Stewartson (DS) equation. The former 
takes the form 

dx [dtu + QudxU + e^dxxxu) + XdyyU = 0, A = ±1 (1.1) 

where (x, y, t) e M^; x Mj^ x Et and where e ^ 1 is a small scaling parameter. The limit 
e — )■ is the dispersionless limit. The case A = — 1 corresponds to the KP I model with 
d. focusing effect, and the case A = 1 corresponds to the KP II model with a defocusing 
effect. The former is also known as the unstable KP equation since the soliton solution 
of the KdV equation is lineary unstable for KP I, whereas it is lineary stable for the 
latter, which is therefore also known as the stable KP equation. These stability 
issues were numerically studied in [?]. It is an interesting result of the present paper 
that both KP equations have very similar numerical convergence properties despite 
completely different stability properties of their exact solutions in [?] . These equations 
appear in different fields of physics in the study of essentially one-dimensional wave 
phenomena with weak transverse effects, for example to model nonlinear dispersive 
waves on the surface of fluids [?]. In this case, KP I is used when the surface tension 
is strong, and KP II when the surface tension is weak. They also model sound 
waves in ferromagnetic media [?], and nonlinear matter- wave pulses in Bose-Einstein 
condensates [?]. The KP equation was introduced by Kadomtsev and Petviashvili in 
[?] to study the stability of the KdV soliton against weak transverse perturbations. 
It was shown to be completely integrable in [?]. Higher dimensional generalizations 
of the KP equations, where the derivative dyy is replaced by the Laplacian in the 
transverse coordinates. A-*- = dyy -\- d^z, are important for instance in acoustics. The 
numerical problems to be expected there are the same as in the 2 -I- 1-dimensional 
case studied here. 

The Davey-Stewartson system can be written in the form 

i£dtu + e^dxxU-ae'^dyyU + 2p(^^+\u\^^u = 0, 

dxx<^> + l3dyy<i> + 2dxx\u\^ = 0, 



where a, /3 and p take the values ±1, where e ^ 1 is again a small dispersion pa- 
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rameter, and where $ is a mean field. Since the e has the same role as the h in the 
Schrodinger equation, the limit e — )■ is also called the semiclassical limit in this 
context. The DS equations are classified [?] according to the ellipticity or hyperbol- 
icity of the operators in the first and second line. The case a = /3 is completely 
integrable [?] and thus provides a 2 + 1-dimensional generalization of the integrable 
NLS equation in 1 + 1 dimensions. The integrable cases are elliptic-hyperbolic called 
DS I, and the hyperbolic-elliptic called DS II. For both there is a focusing {p = —1) 
and a defocusing {p — 1) version. In the following, we will only consider the case DS 
II {a = 1) since the mean field $ is then obtained by inverting an elliptic operator. 
These DS systems model the evolution of weakly nonlinear water waves that travel 
predominantly in one direction, but in which the wave amplitude is modulated slowly 
in two horizontal directions [?], [?]. They are also used in plasma physics [?, ?], to 
describe the evolution of a plasma under the action of a magnetic field. 

Since both KP and DS are completely integrable, there exist many explicit so- 
lutions, which thus provide popular test cases for numerical algorithms. But as we 
will show for the example of KP, these exact solutions, typically solitons, often test 
the equation in a regime where stiffness is not important. The main challenge in 
the study of critical phenomena as dispersive shocks and blowups is, however, the 
numerical resolution of strong gradients in the presence of which the above equations 
are very stiff. This implies that algorithms that perform well for solitons might not 
be efficient in the context studied here. 

Since critical phenomena are generally believed to be independent of the chosen 
boundary conditions, we study a periodic setting Such settings also include rapidly 
decreasing functions which can be periodically continued as smooth functions within 
the finite numerical precision. This allows one to approximate the spatial dependence 
via truncated Fourier series which leads for the studied equations to large stiff systems 
of ODEs, see below. The use of Fourier methods not only gives spectral accuracy in 
the spatial coordinates, but also minimizes the introduction of numerical dissipation 



which is important in the study of dispersive effects. In Fourier space, equations ( 1.1 1 



and (1.2) have the form 

Vt = Lv + N{v,t), (1.3) 

where v denotes the (discrete) Fourier transform of u, and where L and N denote linear 
and nonlinear operators, respectively. The resulting systems of ODEs are classical 
examples of stiff equations where the stiffness is related to the linear part L (it is a 
consequence of the distribution of the eigenvalues of L), whereas the nonlinear part 
contains only low order derivatives. In the small dispersion limit, this stiffness is still 
present despite the small term in L. This is due to the fact that the smaller e 
is, the higher wavenumbers are needed to resolve the rapid oscillations. The first 
numerical studies of exact solutions to the KP equations were performed in [?] and 
[?] and references therein. For the DS system similar studies were done in [?]. In [?] 
blowup for DS was studied for the analytically known blowup solution by Ozawa [?]. 



There are several approaches to deal efficiently with equations of the form (1.31 
with a linear stiff part, implicit-explicit (IMEX), time splitting, integrating factor 
(IF), and deferred correction schemes as well as sliders and exponential time differ- 
encing. To avoid as much as possible a pollution of the Fourier coefficients by errors 
due to the finite difference schemes for the time integration and to allow the use of 



^The boundary conditions will, however, in general influence convergence of the used numerical 
schemes. The restriction in this paper to periodic conditions is due to the studied problems. 
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larger time steps, we mainly consider fourth order schemes. While standard explicit 
schemes impose prohibitively small time steps due to stability requirements (for the 
studied examples the standard fourth order Runge-Kutta (RK) scheme did not con- 
verge for the used time steps), stable implicit schemes are in general computationally 
too expensive in 2 + 1 dimensions. As an example of the latter we consider an implicit 
fourth order Runge-Kutta scheme. The focus of this paper is, however, to compare 
the performance of several explicit fourth order schemes mainly related to exponen- 
tial integrators for various examples in a similar way as in the work by Kassam and 
Trefethen [?] and in [?] for KdV and NLS. 

The paper is organized as follows: In section 2 we briefly list the used numerical 
schemes, integrating factor methods, exponential time differencing, lineary implicit 
schemes, time splitting methods and implicit Runge-Kutta schemes. In section 3 we 
review some analytical facts for the KP equations and study for each of the equations 
an exact solution and an example in the small dispersion limit. In section 4 a similar 
analysis is presented for the semiclassical limit of the focusing and the defocusing DS 
II equation. The found numerical errors are compared in section 5 with the error 
indicated by a violation of the conservation of the L2 norm by the numerical solution. 
In section 6 we add some concluding remarks and outline further directions of research. 

2. Numerical Methods. In this paper we are mainly interested in the numeri- 
cal study of the KP and the DS II equations for Schwartzian initial data in the small 
dispersion limit. The latter implies that we can treat the problem as essentially peri- 
odic, and that we can use Fourier methods. After spatial discretization we thus face 



a system of ODEs of the form (1.3 1. Since we need to resolve high wavenumbers, 



these systems will be in general rather large. The PDEs studied here have high order 



derivatives in the linear part L of (1.3 1, whereas the nonlinear part N contains only 
first derivatives. This means that the stiffness in these systems is due to the linear 
part. The latter will thus be treated with adapted methods detailed below, whereas 
standard methods can be used for the nonlinear part. We restrict the analysis to 
moderate values of the dispersion parameter to be able to study the dependence of 
the different schemes on the time step in finite CPU time. For smaller values of e see 
for instance [?!. 



We will compare several numerical schemes for equations of the form (1.3 1 as in 
[?] and [?] . The PDEs are studied for a periodic setting with periods 2'kLx and 27: Ly 
in X and y respectively. We give the numerical error in dependence of the number 
Nt of time steps as well as the actual CPU time as measured by MATLAB (all 
computations are done on a machine with Intel 'Nehalem' processors with 2.93 GHz 
with codes in MATLAB 7.10). The goal is to provide some indication on the actual 
performance of the codes in practical applications. Since MATLAB is using in general 
a mixture of interpreted and precompiled embedded code, a comparison of computing 
times is not unproblematic. However, it can be done in the present context since 
the main computational cost is due to two-dimensional fast Fourier transformations 
(FFT). For the KP equations all considered approaches (with the exception of the 
Hochbruck-Ostermann ETD scheme which uses 8 FFT commands per time step) use 
6 (embedded) FFT commands per time step as was already pointed out in [?]. For 
the DS II equation, these numbers are doubled since the computation of the mean 
field $ takes another FFT/IFFT pair per intermediate step. Note that an additional 
FFT/IFFT pair per time step is needed in both cases to switch between physical 
and Fourier space (we are interested in a solution in physical space, but the schemes 
are formulated for the Fourier transforms). The (/)-functions in the ETD schemes are 
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also computed via FFT. It can be seen that this can be done with machine precision 
in a very efficient way. Since the (/>-functions have to be obtained only once in the 
computation and since the studied problems are computationally demanding, this only 
has a negligible effect on the total CPU time in the experiments. The numerical error 
is the L2 norm of the difference of the numerical solution and an exact or reference 
solution, normalized by the L2 norm of the initial data. It is denoted by A2. 

2.1. Integrating Factor Methods (IF). These methods appeared first in the 
work of Lawson [?], see [?] for a review. He suggested to take care of the stiff linear 



part of equation (1.3) by using a change of the dep ende nt variables (also called the 
Lawson transformation) w{t) = e~^^v{t). Equation (1.3) becomes 



w'{t) = e-^'N{v,t) (2.1) 

for which we use a fourth order Runge-Kutta (RK) scheme. Hochbruck and Oster- 
mann [?] showed that this IF method has classical order four, but not what they call 
stiff order four. Loosely speaking there can be additional contributions to the error 
of much lower order in the time step due to the stiffness, i.e., in the case of large norm 
||L|| in the integrating factor. They could show that the scheme used here can reduce 
to first order for semilinear parabolic problems. The Hochbruck-Ostermann approach 
uses semigroups. An extension of their theory to hyperbolic equations is possible via 
Co semigroups, which can lead, however, to a slightly lower order of convergence. But 
the results in [?], [?] indicate that similar convergence rates are to be expected for 
hyperbolic equations (both KP and DS are hyperbolic in the sense that the matrix L 



appearing in (1.3 1 has purely imaginary eigenvalues). 



2.2. Driscoll's composite Runge-Kutta Method. The idea of IMEX meth- 
ods (see e.g. [?] for KdV) is the use of a stable implicit method for the linear part of the 



equation (1.3) and an explicit scheme for the nonlinear part which is assumed to be 
non-stiff. In [?] such schemes did not perform satisfactorily for dispersive PDEs which 
is why we only consider a more sophisticated variant here. Fornberg and DriscoU [?] 
provided an interesting generalization of IMEX by splitting also the linear part of 
the equation in Fourier space into regimes of high, medium, and low wavenumbers, 
and by using adapted numerical schemes in each of them. They considered the NLS 
equation as an example. Driscoll's [?] idea was to split the linear part of the equation 
in Fourier space just into regimes of high and low wavenumbers. He used the fourth 
order RK integrator for the low wavenumbers and the lineary implicit RK method 
of order three for the high wavenumbers. He showed that this method is in practice 
of fourth order over a wide range of step sizes. We confirm this here for the cases 
where the method converges, which it fails to do, however, sometimes for very stiff 
problems. In particular, he used this method for the KP II equation at the two phase 
solution we will also discuss in this paper as a test case. We call the method DCRK 
in the following. 

2.3. Exponential Time Differencing Methods. Exponential time differenc- 
ing schemes were developed originally by Certaine [?] in the 60s, see [?] and [?] for 
comprehensive reviews of ETD methods and their history. The basic idea is to use 



equidistant time steps h and to integrate equation (1.3) exactly between the time 



steps tn and tn+i with respect to t. With v{tn) — v„ and v{t„+i) ~ Vn+i, we get 

'vn + / e^(''-^)N(w(t„ + r), t„ + T)dT. 
Jo 



The integral will be computed in an approximate way for which different schemes 
exist. Wc use here only Runge-Kutta schemes of classical order 4, Cox-Matthews [?], 
Krogstad [?] and Hochbruck-Ostermann [?]. The latter showed that the stiff order 
of the Cox-Matthews scheme is only two, and the one of Krogstad's is three. Notice 
that both schemes can be of stiff order four if the studied system satisfies a certain 
number of non-trivial auxiliary conditions, see [?]. As the numerical tests show in 
the following, order reduction can be observed in some cases, but not in all. We will 
speak in the following of the stiff regime of an equation where order reduction can 
be observed, or where certain schemes do not converge, and of the non-stiff regime, 
where this is not the case. Note that this is not a rigorous definition. Both Cox- 
Matthews' and Krogstad's schemes are, however, four-stage methods, whereas the 
Hochbruck-Ostermann method is a five-stage method that has stiff order four. Thus 
all these methods should show the same convergence rate in the non-stiff regime, but 
could differ for some problems in the stiff regime. Notice that these results [?] were 
established for parabolic PDEs, and that the applicability for hyperbolic PDEs of the 
type studied here is not obvious. One of the purposes of our study is to get some 
experimental insight whether the Hochbruck-Ostermann theory holds also in this case. 

The main technical problem in the use of ETD schemes is the efficient and accurate 
numerical evaluation of the functions 



i.e., functions of the form (e^ — l)/-^ and higher order generalizations thereof, where 
one has to avoid cancellation errors. Kassam and Trefethen [?] used complex contour 
integrals to compute these functions. The approach is straight forward for diagonal 
operators L that occur here because of the use of Fourier methods: one considers a unit 
circle around each point z and computes the contour integral with the trapezoidal rule 
which is known to be a spectral method in this case. Schmelzer [?] made this approach 
more efficient by using the complex contour approach only for values of z close to the 
pole, e.g with \z\ < 1/2. For the same values of z the functions 0^ can be computed 
via a Taylor series. These two independent and very efficient approaches allow a 
control of the accuracy. We find that just 16 Fourier modes in the computation of the 
complex contoiir integral are sufficient to determine the functions (jji to the order of 
machine precision. Thus we avoid problems reported in [?], where machine precision 
could not be reached by ETD schemes due to inaccuracies in the determination of 
the (^-functions. The computation of these functions takes only negligible time for 
the 2 -|- 1-dimensional equations studied here, especially since it has to be done only 
once during the time evolution. We find that ETD as implemented in this way has 
the same computational costs as the other used schemes. 

2.4. Splitting Methods. Splitting methods are convenient if an equation can 
be split into two or more equations which can be directly integrated. The motivation 
for these methods is the Trotter-Kato formula [?, ?] 



where A and B are certain unbounded linear operators, for details see [?]. In particular 
this includes the cases studied by Bagrinovskii and Godunov in [?] and by Strang [?]. 
For hyperbolic equations, first references are Tappert [?] and Hardin and Tappert [?] 
who introduced the split step method for the NLS equation. 





(2.2) 
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The idea of these methods for an equation of the form ut = {A + B) u is to write 
the solution in the form 

u{t) — exp(cii^) exp{ditB) exp(c2tA) exp{d2tB) ■ ■ ■ exp{cf^tA) exp{dktB)u{0) 

where (ci, . . . , Ck) and {di, . . . , dk) are sets of real numbers that represent fractional 
time steps. Yoshida [?] gave an approach which produces split step methods of any 
even order. 

The KP equation can be split into 

Ut + 6uux — 0, (2-3) 
(J"M)( - iklJ'iu] + A— ^ -FM = 0, (2.4) 

where here and in the following we write the 2-dimensional Fourier transform of u in 
the form 

J'iu]:^ I u{x,y,t)e-'''''''~'^yydxdy. (2.5) 



The Hopf equation (2.3) can be integrated in implicit form with the method of char- 



acteristics, and the linear equation in Fourier space (2.4 1 can be directly integrated, 
but the implicit form of the solution of the former makes an iteration with interpola- 
tion to the characteristic coordinates necessary that is computationally too expensive. 
Therefore we consider splitting here only for the DS equation. The latter can be split 
into 

ieut = e^{-Uxx + auyy), ^xx + ct^yy + 2 (\u\^) =0, (2.6) 

V / XX 

-2p[^ + \uf^u, (2.7) 



leut 



which are explicitly integrable, the first two in Fourier space, equation (2.7) in physical 
space since |up is a constant in time for this equation. Convergence of second order 
splitting along these lines was studied in [?]. We study here second and fourth order 
splitting schemes for DS as given in [?]. 

2.5. Implicit Runge Kutta Scheme. The general formulation of an s-stage 
Runge-Kutta method for the initial value problem y' = f{y,t), y(to) — yo is the 
following: 



yn+l =yn + biKi, (2.8) 



K, = f[tn + c^h, yn + /i^ aijKj ) , (2.9) 

s 

where bi, Uij, i, j = 1, s are real numbers and q = ^ aij. 

For the implicit Runge-Kutta scheme of order 4 (IRK4) used here (Hammer- 

HoUingsworth method), we have ci = |— C2 = ^ + '^ii — '^22 = 1/4, 
ai2 = \ — ^21 = i + ^ s-iid &i = &2 = 1/2. This scheme is also known as the 
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2-stage Gauss method. It is of classical order 4, but stage order 2. This is the reason 
why an order reduction to second order can be observed in certain examples. 

The implicit character of this method requires the iterative solution of a high 
dimensional system at every step which is done via a simplified Newton method. For 



the studied examples in the form (1.3), we have to solve equations of the form 



y^Ay + b{y) 

for y, where A is a linear operator independent of y, and where 6 is a vector with a 
nonlinear dependence on y. These are solved iteratively in the form 

Vn+i = (1 - A)-i&(2/„). 

By treating the linear part that is responsible for the stiffness explicitly as in an IMEX 
scheme, the iteration converges in general quickly. Without taking explicit care of the 
linear part, convergence will be extremely slow. The iteration is stopped once the Loo 
norm of the difference between consecutive iterates is smaller than some threshold 
(in practice we work with a threshold of 10~^). Per iteration the computational cost 
is essentially 2 FFT/IFFT pairs. Thus the IRK4 scheme can be competitive with 
the above explicit methods which take 3 or 4 FFT/IFFT pairs per time step if not 
more than 2-3 iterations are needed per time step. This can happen in the below 
examples in the non-stiff regime, but is not the case in the stiff regime. We only test 
this scheme where its inclusion appears interesting and where it is computationally 
not too expensive. 

3. Kadomtsev-Petviashvili Equation. In this section we study the efficiency 
of the above mentioned numerical schemes in solving Cauchy problems for the KP 
equations. We first review some analytic facts about KP I and KP II which are impor- 
tant in this context. Since the KP equations are completely integrable, exact solutions 
exist that can be used as test cases for the codes. We compare the performance of 
the codes for the exact solutions and a typical example in the small dispersion limit. 

3.1. Analytic Properties of the KP Equations. We will collect here some 
analytic aspects of the KP equations which will be important for an understanding 
of several issues in the numerical solution of Cauchy problems for the KP equations, 
see [?] for a recent review and references therein. 

In this paper we will look for KP solutions that are periodic in x and y, i.e., for 
solutions on xK. This includes for numerical purposes the case of rapidly decreasing 
functions in the Schwartz space iS(M^) if the periods are chosen large enough that |m| is 
smaller than machine precision (we work with double precision throughout the article) 
at the boundaries of the computational domain. Notice, however, that solutions to 
Cauchy problems with Schwartzian initial data UQ{x,y) will not stay in 5(M^) unless 
UQ(x,y) satisfies an infinite number of constraints. This behaviour can be already 
seen on the level of the linearized KP equation, see e.g. [?, ?], where the Green's 
function implies a slow algebraic decrease in y towards infinity. This leads to the 
formation of tails with an algebraic decrease to infinity for generic Schwartzian initial 
data. The amplitude of these effects grows with time (see for instance [?]). In our 
periodic setting this will give rise to echoes and a weak Gibbs phenomenon at the 
boundaries of the computational domain. The latter implies that we cannot easily 
reach machine precision as in the KdV case unless we use considerably larger domains. 
As can be seen from computations in the small dispersion limit below and the Fourier 
coefficients in sect. 5, we can nonetheless reach an accuracy of better than lO^^'^ on 



the chosen domain. For higher precisions and larger values of t, the Gibbs phenomena 
due to the algebraic tails become important. 

The KP equation ([lT]) is not in the standard form for a Cauchy problem. As 
discussed in [?] t is not a tinielike but characteristic coordinate if the dispersionless 
KP equation (e = 0) is considered as a standard second order PDE. In practice one 



is, however, interested in the Cauchy problem for t. To this end one writes (1.1) in 
evolutionary form 

dfU + 6udxU + e^dxxxU — —Xd^^dyyU, A = ±1. (3-1) 



Equations (3.1) and (1.1) are equivalent for certain classes of boundary conditions 
as periodic or rapidly decreasing at infinity. Since we will always impose periodic 
boundary conditions in the following, both forms of the KP equation are equivalent 
for our purposes. The antiderivative is to be understood as the Fourier multiplier 
with the singular symbol —i/kx- In the numerical computation this multiplier is 
regularized in standard way (similar to the Dirac regularization of 1/x by adding an 
arbitrary small imaginary part iO to x) as 



kx + iXd ' 

where we choose S = eps = 2~^^ ~ 2.2 * 10^^^. Since the typical precision to be 
achieved in the studied examples with double precision is of the order 10^^** because 
of rounding errors, this is essentially equivalent to adding a numerical zero (see also 
the discussion in [?]). 

The divergence structure of the KP equations has the consequence that 

dyyu{x,y,t)dx ^0, Vi > 0, (3.2) 

T 

even if this constraint is not satisfied for the initial data uo{x,y). It was shown in 
[?, ?] that the solution to a Cauchy problem not satisfying the constraint will not be 
smooth in time for t = 0. Numerical experiments in [?] indicate that the solution 
after an arbitrary small time step will develop an infinite 'trench' the integral over 



which just ensures that (3.2) is fulfilled. To propagate such initial data, the above 
regularization is in fact needed (for data satisfying the constraint via the condition 
that the Fourier coefficients for kx = vanish, this property could be just imposed at 
each time step). The infinite trench due to initial data not satisfying the constraint 
implies a rather strong Gibbs phenomenon. To avoid the related problems, we always 



consider initial data that satisfy (3.2). A possible way to achieve this is to consider 
data that are x-derivatives of periodic or Schwartzian functions. 

The complete integrability of the KP equations implies that efficient tools exist 
for the generation of exact solutions. We always put e = 1 for the exact solutions. For 
a recent review of the integrable aspects of KP see [?] . The most popular KP solutions 
are line solitons, i.e., solutions localized in one spatial direction and infinitely extended 
in another. Such solutions typically have an angle not equal to or 90 degrees with 
the boundaries of the computational domain, which leads to strong Gibbs phenomena. 
This implies that these solutions are not a good test case for a periodic setting. If 
the angle is or 90 degrees, the solution only depends on one of the spatial variables 
and thus does not test a true 2d code. There exists a lump soliton for KP I which is 
localized in all spatial directions, but only with algebraic fall off. This would again 
lead to strong Gibbs phenomena in our setting. 
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However a solution due to Zaitsev [?] to the KP I equation is localized in one 
direction and periodic in the second (a transformation of the form x ^ ix, y ^ iy 
exchanges these two directions). It has the form 



u{i,y)^2a- — — ——2 (3.3) 



(cosh(a^) — (3cos{Sy)) 

where 



(^x-ct, c^a ^— ^, and S ^ ^ ^—^a 

This solution is localized in x, periodic in y, and unstable as discussed in [?]. 

Algebro-geometric solutions to the KP equation can be constructed on an arbi- 
trary compact Riemann surface, see e.g. [?], [?]. These solutions are in general almost 
periodic. Solutions on genus 2 surfaces, which are all hyperelliptic, are exactly peri- 
odic, but in general not in both x and y. A doubly periodic solution with exactly this 
property of KP II of genus 2 can be written as 

(92 

uix,y,t) = 2—lii9{ip,,ipr,B) (3.4) 
where 9 {ipi^ ip2] B) is defined by the double Fourier series 

oo oo 

e{^i,^2;B)^ J2 e5™^^™+""^^ (3.5) 

mi— — oom2— — oo 

where rn^ = (mi, m2), and where _B is a 2 x 2 symmetric, negative-definite Riemann 
matrix 

_B = ^ 6A^'+ d ) ' ^^^^ ^^^^ parameters X ^ 0, b and d. 

The phase variable ip has the form tpj = fijX + Vjy + ujjt + ipjfi, j = 1,2. The solution 
travels as the Zaitsev solution with constant speed in x-direction. 

Remark 3.1. The standard 4th order Runge-Kutta scheme did not converge for 
any of the studied examples for the used time steps. The reason is that the Fourier 
multiplier —i/k^ imposes very strong stability restrictions on the scheme. 

3.2. Numerical solution of Cauchy problems for the KP I equation. 



Zaitsev solution. We first study the case of the Zaitsev solution (3.3) with a = 1 
and f3 — 0.5. Notice that this solution is unstable against small perturbations as 
shown numerically in [?], but that it can be propagated with the expected numerical 
precision by the used codes. As initial data we take the solution centered at —Lx/2 
(we use — Ly = 5) and propagate it until it reaches The computation is 

carried out with 2^^ x 2^ points for {x^y) S [— Stt, 57r] x [— Stt, Stt] and t < 1. The 



decrease of the numerical error is shown in Fig. 3.1 in dependence of the time step and 



in dependence of the CPU time. A linear regression analysis in a double logarithmic 



plot (logiQ A2 = — alogiQ Nt + b) is presented in Fig. 3.1 where we can see that all 
schemes show a fourth order behavior: we find a = 4.32 for the Integrating Factor 
method, a = 4.38 for DCRK method, a = 3.93 for Krogstad's ETD scheme, a = 4 
for the Cox-Matthews scheme, and a = 3.98 for the Hochbruck-Ostermann scheme. 
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Fig. 3.1. Normalized L2 norm of the numerical error in the solution to the KP I equation with 
initial data given by the Zaitsev solution for several numerical methods, as a function of Nt (left) 
and as a function of CPU time (right). 



In this context the DCRK method performs best, foUowed by the ETD schemes that 
have almost identical performance (though the Hochbruck-Ostermann method uses 



more internal stages and thus more CPU time in Fig. 3.1 1. It can also be seen that 
the various schemes do not show the phenomenon of order reduction as discussed in 
[?] , which implies that the Zaitsev solution tests the codes in a non-stiff regime of the 
KP I equation. 

Small dispersion limit for KP I. To study KP solutions in the limit of small 
dispersion (e — ^ 0), we consider Schwartzian initial data satisfying the constraint 



(3.2). As in [?] we consider data of the form 



^0(2:, y) = — 9a;Sech^(i?) where R = \J + , (3,g) 



By numerically solving the dispersionless KP equation (put e — in (1.1)), we de- 
termine the critical time of the appearance of a gradient catastrophe by the breaking 
of the code, see [?]. To study dispersive shocks, we run the KP codes for some time 



larger than this critical time. The solution can be seen in Fig. 3.2 It develops tails 
with algebraic fall off towards infinity. The wave fronts steepen on both sides of the 
origin. In the regions of strong gradients, rapid modulated oscillations appear. For a 
detailed discussion, see [?]. 

The computation is carried out with 2^^ x 2^ points for (a;,y) G [— Stt, Stt] x 
[— 57r, Stt], e = 0.1 and t < 0.4. As a reference solution, we consider the solution 
calculated with the Hochbruck-Ostermann method with Nt = 5000 time steps. The 
normalized L2 norm of the difference between this reference solution and the numerical 
solution is shown in Fig. |3.3| in dependence on the time step with a regression analysis 
and in dependence on the CPU time. Here we can see clearly the phenomenon of order 
reduction established analytically for parabolic systems by Hochbruck and Ostermann 
[?]. In the stiff regime (here up to errors of order 10"'*) DCRK does not converge, the 
Integrating Factor method shows only first order behaviour (as predicted in [?]), the 
IRK4 scheme shows second order convergence, and ETD methods perform best. This 
implies that the Cox-Matthews and Krogstad method with similar performance are 
the most economic for the stiff regime of the KP I equation, which gives the precision 
one is typically interested in in this context. For higher precisions we find a = 1.87 
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t=0.0 



t=0.16 





t=0.4 




Fig. 3.2. Solution to the KP I equation for the initial data uq = —dxsech^ (R) where R 



\/ + for several values of t. 
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Fig. 3.3. Normalized L2 norm of the numerical error for the solution shown in Fig. \3.2\ for 
several numerical methods, as a function of Nt (left) and as a function of CPU time (right). 



for the Integrating Factor method, a = 4.20 for DCRK, a = 4.03 for IRK4, a = 3.75 
for Krogstad's ETD scheme, a = 3.90 for the Cox-Matthews scheme, and a — 3.96 for 
the Hochbruck-Ostermann scheme. 

To study empiricahy the phenomenon of order reduction in exponential integra- 
tors, and to observe the transition from a stiff to a non stiff regime we study the ETD 



schemes in more detail in Fig. 3.4 This is indicated by the fact that ETD schemes are 

12 



only of order three in this stiff region instead of order four. It appears that all schemes 




Fig. 3.4. Phenomenon of order reduction in exponential integrators for KP I and KP II re- 
spectively in the small dispersion limit, ETD schemes of Fig. \3. 3\ (left) and of Fig. \3. 7| ( right). 

show a slight order reduction though this is not the case for the Hochbruck-Osterniann 
method in the parabolic case. 

3.3. Numerical solution of Cauchy problems for the KP II equation. 

Doubly periodic .solution of KP II. The computation for the doubly periodic so- 
lution to KP II is carried out with 2* x 2^ points for {x,y) S [— 57r, Stt] x [— Stt, Stt] 
and t < 1 with the parameters 6 = 1, A = 0.15, bX^ + d = —1, /ii = /i2 = 0.25, 
= -ly^ = 0.25269207053125, uji = uj2 ^ -1.5429032317052, and tpi^o = V2,o = 0. 
The decrease of the numerical error is shown in Fig. |3.5| in dependence of Nt and in 
dependence on CPU time. From a linear regression analysis in a double logarithmic 




Fig. 3.5. Normalized L2 norm of the numerical error for the time evolution of the doubly 
periodic solution to the KP II equation for several numerical methods, as a function of Nt (left) and 
as a function of GPU time (right). 

plot we can see that all schemes are fourth order: one finds a = 4.45 for the Integrating 
Factor method, a = 4.33 for DCRK, a = 3.97 for Krogstad's ETD scheme, a = 4 for 
the Cox-Matthews scheme, and a — 3.95 for the Hochbruck-Ostermann scheme. As 
for the Zaitsev solution, DCRK performs best followed by the ETD schemes. We thus 
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confirm DriscoU's results in [?] on tlie efficiency of his metliod for this example. The 
absence of order reductions indicates again that the exact solution tests the equation 
in a non-stiff regime. IRK4 is competitive for larger time steps in this case since only 
very few iterations (1-3) are needed. 

Small dispersion limit for KP II. We consider the same initial data and the same 
methods as for KP I. In Fig. |3.6| the time evolution of these data can be seen. The 
solution develops tails this time in negative a;-direction. The steepening of the wave 
fronts happens at essentially the same time, but the gradients are stronger in the 
vicinity of the tails (see [?]). This is also where the stronger oscillations appear. 



t=0.0 



t=0.16 




t=0.18 





Fig. 3.6. Solution to the KP 11 equation for the initial data uq = —dxSech?{R) where R = 
X'^ + y'^ for several values oft. 

The computation is carried out with 2^^ x 2^ points for {x,y) G [—5tt, Stt] x 
[— 57r, Stt], e = 0.1 and t < 0.4. As a reference solution, we consider the solution 
calculated with the Hochbruck-Ostermann method with Nt = 5000 time steps. The 
dependence of the normalized L2 norm of the difference between this reference solution 
and the numerical solution on Nt and on CPU time is shown in Fig. |3.7[ We obtain 
similar results as in the small dispersion limit of KP I: for typical accuracies one 
is interested in in this context, DCRK does not converge, the Integrating Factor 
method shows only a first order behavior, the IRK4 method is of second order, and 
the ETD methods perform best. In the non-stiff regime we find a = 1.06 for the 
Integrating Factor method, a = 4.42 for DCRK, a = 4.33 for IRK4, a = 4.15 for 
Krogstad's ETD scheme, a = 4.07 for the Cox-Matthews scheme, and a — 3.99 for 
the Hochbruck-Ostermann scheme. Once again, we study empirically the phenomenon 
of order reduction in exponential integrators, and observe a transition from a stiff to 
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a non stiff region (Fig. 3.3 ), indicated by the fact that ETD schemes are only of order 



three in this stiff region instead of of order four. 

4. Davey-Stewartson equation. In this section we perform a similar study as 
for KP of the efficiency of fourth methods in solving Cauchy problems for the DS II 
equations. We first review some analytic facts about the focusing and defocusing DS 
II equations which are of importance in this context. We compare the performance 
of the codes for a typical example in the small dispersion limit. 

4.1. Analytic properties of the DS equations. For a review see for instance 
the book by Sulem and Sulem [?]. We will study here only the DS II equations (a = 



P — 1 in eq. ( 1.2 )) since the elliptic operator for $ can be inverted by imposing simple 
boundary conditions. For a hyperbolic operator acting on $, boundary conditions for 
wave equations have to be used. 



We will consider the equations again on T 



Due to the ellipticity of the 



operator in the equation for it can be inverted in Fourier space in standard manner 
by imposing periodic boundary conditions on $ as well. As before this case contains 
Schwartzian functions that are periodic for numerical purposes. Notice that solutions 
to the DS equations for Schwartzian initial data stay in this space at least for finite 
time in contrast to the KP case. Using Fourier transformations $ can be eliminated 



from the first equation by a transformation of the second equation in (1.2) and an 
inverse transformation. With (2.5) we have 



-2J- 



1.2 



Zr-2 



which leads in (1.2) as for KP to a nonlocal equation with a Fourier multiplier. This 



implies that the DS equation requires an additional computational cost of 2 two- 
dimensional FFT per intermediate time step, thus doubling the cost with respect to 
the standard 2d NLS equation. Notice that from a numerical point of view the same 
applies to the elliptic-elliptic DS equation that is not integrable. Our experiments 
indicate that except for the additional FFT mentioned above, the numerical treatment 
of the 2d and higher dimensional NLS is analogous to the DS II case studied here. The 
restriction to this case is entirely due to the fact that one can hope for an asymptotic 
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description of the small dispersion limit in the integrable case. Thus we study initial 
data of the form UQ{x,y) = a{x,y)exp{ib{x,y)/e) with a, 6 e M, i.e., the semi-classical 
limit well known from the Schrodinger equation. Here we discuss only real initial data 
for convenience. 

It is known that DS solutions can have blowup. Results by Sung [?] establish 
global existence in time for initial data ^po & Lp, I < p < 2 with a Fourier transform 
J-lipo] S Li n Loo subject to the smallness condition 

mM\LA\^m\\L^ < Y (^^) ^^-^^ 

in the focusing case. There is no such condition in the defocusing case. Notice that 



condition (4.1 ) has been established for the DS II equation with e = 1. The coordinate 



change x' — x/e, t' — t/e transforms the DS equation (1.2) to this standard form. 
This implies for the initial data uq = exp{—x'^ — r]y^) we study for the small dispersion 



limit of the focusing DS II system in this paper that condition (4.1) takes the form 

0.0477. 




This condition is not satisfied for the values of e and rj we use here. Nonetheless we 
do not observe any indication of blowup on the shown timescales. One of the reasons 
is that the rescaling with e above also rescales the critical time for blowup by a factor 
1/e. In addition it is expected that the dispersionless equations will for generic initial 
data have a gradient catastrophe at some time < oo, and that the dispersion will 
regularize the solution for small times t > tc > 0. However there are no analytic 
results in this context. 

The complete integrability of the DS II equation implies again the existence of 
explicit solutions. Multi-soliton solutions will be as in the KP case localized in one 
spatial direction and infinitely extended in another, the lump solution is localized in 
two spatial directions, but with an algebraic fall off towards infinity. Thus these are 
again not convenient to test codes based on Fourier methods as in the KP case. Since 
the study of the small dispersion limit below indicates that the time steps have to be 
chosen sufficiently small for accuracy reasons such that in contrast to KP no order 
reduction observed, we will not study any exact solutions here. 

4.2. Small dispersion limit for DS II in the defocusing case. We consider 
initial data uq of the form 

uo(x, y) — e^^ , where R = + 77?;^ with = 1 (4-2) 

and use the same methods as before together with time splitting methods of order 2 
and one of order 4, as explained in section 2.4. The defocusing effect of the defocusing 



DS II equation for these initial data can be seen in Fig. 4.1 where jwp is shown for 
several values of t. The compression of the initial pulse into some almost pyramidal 
shape leads to a steepening on the 4 sides parallel to the coordinate axes and to 
oscillations in these regions. 



The computations are carried out with 2 x 2 " points for {x,y) G [— Stt, Stt] 



X 



[— Stt, Stt], e — 0.1 and t < 0.8. To determine a reference solution, we compute 
solutions with 6000 time steps with the ETD, the DCRK and the IF schemes and take 
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t=0.28 




Fig. 4.1. Solution to the defocusing DS II equation for the initial data no = exp(—R^) where 
R = \J -\-'ip- and t = 0.1 for several values of t. 

the arithmetic mean. The dependence of the normahzed L2 norm of the difference of 
the numerical sohitions with respect to this reference solution on Nf and on CPU time 
is shown in Fig. |4.2[ A linear regression shows that all fourth order schemes show a 




Fig. 4.2. Normalized L2 norm of the numerical error for several numerical methods for the 
situation shown in Fig. \4 l\ as a function of Nt (left) and of CPU time (right). 

fourth order behavior except for IRK4 (a — 1.78), as is obvious from the straight lines 
with slope a = 3.95 for the Integrating Factor method, a = 4.01 for DCRK, a — 3.99 
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for Krogstad's ETD scheme, a — 3.99 for the Cox-Matthews scheme, a = 3.95 for 
the Hochbruck-Ostermann scheme, and a = 3.86 for the time spUtting method. The 
second order spUtting scheme shows the expected convergence rate and performs very 
weh for lower precision. For smaller time steps, the advantage of the fourth order 
schemes is more pronounced. Apparently the system is 'stiff' for the IRK4 scheme 
since it only shows second order behavior. Notice that the time splitting scheme 
reaches its maximal precision around 10^^, a behavior which was already noticed in 
[?] for the study of the Nonlinear Schrodinger equation in the small dispersion limit. 
It appears that this behavior is due to resonances of errors of the split equations, but 
the identification of the precise reason will be the subject of further research. The 
same effect is observed for second order splitting for smaller time steps than shown 
in Fig. |4.2[ However, both schemes work very well at the precisions in which one is 
normally interested in. We mainly include a second order scheme here because of the 
additional computational cost due to the function <i> in the DS system. This could 
make a second order scheme competitive in terms of CPU time because of the lower 



number of FFT used per time step. It can be seen in Fig. 4.2 that this is not the case. 



We conclude that the ETD schemes perform best in this context. 

4.3. Small dispersion limit for the focusing DS II equation. For the 



focusing DS II in the small dispersion limit we consider initial data of the form (4.2 1 
with 77 ~ 0.1 and the same methods as before. The focusing effect of the equation 
can be clearly recognized in Fig. |4.3| The initial peak grows until a breakup into a 
pattern of smaller peaks occurs. 



t=0.0 t=0.276 




Fig. 4.3. Solution to the focusing DS II equation for the initial data uq = 
exp(— i?^) where R = \J xP- + r]'ip- , X] = 0.1 and e = 0.1 for several values oft. 
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It is crucial to provide sufficient spatial resolution for the central peak. As for 
the 1+1-dimensional focusing NLS discussed in [?], the modulational instability of 
the focusing DS II leads to numerical problems if there is no sufficient resolution for 
the maximum. In [?] a resolution of 2^'^ modes was necessary for initial data 
and e = 0.1 for the focusing NLS in 1 + 1 dimensions. The possibility of blowup 
in DS requires at least the same resolution despite some regularizing effect of the 
nonlocality <&. With the computers we could access, a systematic study of time 
integration schemes with a resolution of 2^^ x 2^^ was not possible in Matlab. Thus 
we settled for initial data close to the one-dimensional case, which allowed for a lower 
resolution, see the next section for the Fourier coefficients. The computation is carried 
out with 2^2 X 2" points for e [-Stt, Svr] x [-Stt, Stt], e = 0.1 and t < 0.6. To 

determine a reference solution, we compute solutions with 6000 time steps with the 
ETD, the DCRK and the IF schemes and take the arithmetic mean. The dependence 
of the normalized L2 norm of the difference of the numerical solutions with respect 



to this reference solution on Nt and on CPU time is shown in Fig. 4.4 All schemes 




Fig. 4.4. Normalized L2 norm of the numerical error for the example of Fig. \4-3\ for several 
numerical methods as a function of Nt (left) and of GPU time (right). 

except the time splitting scheme show a fourth order behavior, as can be seen from 
the straight lines with slope a = 3.36 for the Integrating Factor method, a = 3.93 for 
DCRK, a — 4.06 for Krogstad's ETD scheme, a — 4.11 for the Cox-Matthews scheme, 
a = 4.13 for the Hochbruck-Ostermann scheme, and a = 2.5 for the fourth order time 
splitting method. We conclude that in this context DCRK performs best, followed by 
the ETD schemes. We do not present results for the IRK4 scheme here since it was 
computationally too expensive. 

5. Numerical conservation of the L2 norm. The complete integrability of 
the KP and the DS equations implies the existence of many or infinitely many con- 
served quantities (depending on the function spaces for which the solutions are de- 
fined). It can be easily checked that the Li norm and the L2 norm of the solution 
are conserved as well as the energy. We do not use here symplectic integrators that 
take advantage of the Hamiltonian structure of the equations. Such integrators of 
fourth order will be always implicit which will be in general computationally too ex- 
pensive for the studied equations as the experiment with the implicit IRK4 scheme 
showed. Moreover it was shown in [?] that fourth order exponential integrators clearly 
outperform second order symplectic integrators for the NLS equation. 
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The fact that the conservation of L2 norm and energy is not implemented in 
the code ahows to use the 'numerical conservation' of these quantities during the 
computation or the lack thereof to test the quality of the code. We will study in this 
section for the previous examples to which extent this leads to a quantitative indicator 



of numerical errors. Note that due to the non-locality of the studied PDEs ( 1.1 ) and 



(1.2), the energies both for KP, 

E[u[t)] ■.= \ ( {dMt, X, v)f - \[d-^dyu{t, X, y)f - 2e\^{t, x, y)) dxdy, 



and for DS II, 



'\d^u{t,x,y)\^ ~ ^\dyu{t,x,y)'{^ 



-p {\u{t, X, y)\^-\ {Ht. ^, Vf + {d-'dyHt. X, y)f)^ 



dxdy, 



contain anti-derivatives with respect to x. Since the latter are computed with Fourier 
methods, i.e., via division by kx in Fourier space, this computation is in itself numeri- 
cally problematic and could indicate problems not present in the numerical solution of 
the Cauchy problem. Therefore we trace here only the L2 norm \u(t,x,y)\'^dxdy, 
where these problems do not appear. In the plots we show the variable test de- 
fined as test — M{t)/M{{)) — 1, where M{t) is the numerically computed L2 norm in 
dependence of time. 

Notice that numerical conservation of the L2 norm can be only taken as an indi- 
cation of the quality of the numerics if there is sufficient spatial resolution. Therefore 
we will always present the Fourier coefficients for the final time step for the considered 
examples. No dealiasing techniques are used. We will discuss below the results for 
the small dispersion limit. 



For the KP I example of Fig. 3.2 we get the Fourier coefficients at the final 



time and the mass conservation shown in Fig. 5.1 It can be seen that the Fourier 
coefficients decrease in A:a;-direction to almost machine precision, whereas this is not 
fully achieved in fcj,-direction. This is partly due to the necessity to allow extensive 
studies of the dependence on the time-stepping in finite computing time and thus to 
keep the spatial resolution low, and partly due to a Gibbs phenomenon mainly in 



/cy-direction due to the formation of the algebraic tails in Fig. 3.2 Mass conservation 



can be seen to be a viable indicator of the numerical accuracy by comparing with 



Fig. 3.3 in the range of accuracy in which one is typically interested 10"**), mass 
conservation overestimates the actual accuracy by roughly 2 orders of magnitude. It 
can be seen that it shows also at least a fourth order decrease. 



The situation is very similar for the small dispersion example for KP II of Fig. 3.6 
as can be seen in Fig. |5.2[ 

For the defocusing DS II equation and the example shown in Fig. |4.1[ the Fourier 
coefficients decrease to machine precision despite the lower resolution than for KP. One 
reason for this is the absence of algebraic tails in the solution. The mass shows as for 
KP at least fourth order dependence on the time step and overestimates the numerical 
precision by roughly two orders of magnitude. This is not true for the splitting scheme 
for which mass conservation is no indication of the numerical precision at all. This 



seems to be due to the exact integration of the equations (2.7) into which DS is split 
(for one of them the L2 norm is constant). The found numerical mass does not appear 
to reflect the splitting error that is the reason for the numerical error here. 

20 



Fig. 5.1. No n-con servation of the numerically computed L2 norm of the solution to the problem 
considered in Fig. \3.3\ in dependence on the time step (left) and the Fourier coefficients for the final 
time (right). 




Fig. 5.2. Non-conservation of the numerically computed L2 norm of the solution to the problem 
considered in Fig. \3. T\ in dependence on the time step (left) and the Fourier coefficients for the final 
time (right). 



For the small dispersion example for the focusing DS II equation of Fig. |4.1| it can 
be seen in Fig. |5.4| that spatial resolution is almost achieved. There is a certain lack of 
resolution in the direction which leads to the formation of some structure close to 
ky = 0. This is related to the modulational instability of solutions to the focusing DS 
II equation. It will disappear for higher resolutions. Numerical conservation of the 
L2 norm of the solution overestimates numerical accuracy by 2-3 orders of magnitude 
for an error of the order of 10"'^. Once more it cannot be used as an indicator for 
the numerical error in the splitting case, where it is almost independent of the time 
step. For the other cases numerical conservation of the L2 norm shows a dependence 
on Nt between fourth and fifth order. This indicates as for the NLS case in [?] that 
the numerical error has a divergence structure which leads to a higher order decrease 
of the L2 norm than for the actual error. This behavior is also present in the above 
examples, but less pronounced. 
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Fig. 5.3. No n-con servation of the numerically computed L2 norm of the solution to the problem 
considered in Fig. \4-2\ in dependence on the time step (left) and the Fourier coefficients for the final 
time (right). 




Fig. 5.4. Non-conservation of the numerically computed L2 norm of the solution to the problem 
considered in Fig. \4.4\ ™ dependence on the time step (left) and the Fourier coefficients for the final 
time (right). 

6. Conclusion. It was shown in this paper that fourth order time stepping 
schemes can be efficiently used for higher dimensional generalizations of the KdV and 
the NLS equations, where the stiffness of the system of ODEs obtained after spatial 
discretization can be a problem. Implicit schemes as IRK4 are computationally too 
expensive in the stiff regime, whereas standard explicit schemes as RK require for 
stability reasons too restrictive requirements on the time steps for the KP and DS 
equations. For these equations the non-localities in the PDEs lead to singular Fourier 
multipliers which make standard explicit schemes in practice unusable for stability 
reasons. 

IMEX schemes do not converge in general for similar reason. DriscoU's composite 
RK variant is generally very efficient if the studied system is not too stiff, but fails 
to converge for strong stiffness. Exponential integrators do not have this problem. 
The order reduction phenomenon is a considerable problem for IF schemes in the stiff 
regime, but less so for ETD schemes. The Hochbruck-Ostermann method performs 
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in general best, but the additional stage it requires is in practice not worth the effort 
in comparison with Krogstad's or Cox-Matthews' method. The computation of the 
(/)-functions in ETD is inexpensive for the studied problems since it has to be done 
only once. 

Since stiffness is not the limiting factor for DS II, all schemes perform well in this 
context. But the modulational instability of the focusing case requires high spatial 
resolution we could not achieve in Matlab on the used computers for more general 
initial data. Thus the code will be parallelized to allow the use of higher spatial 
resolution without allocating too much memory per processor. 
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